Epicardial Dispersion of Repolarization Promotes the Onset of Reentry in Brugada Syndrome: A Numerical Simulation Study

The Brugada syndrome (BrS) is a cardiac arrhythmic disorder responsible for sudden cardiac death associated with the onset of ventricular arrhythmias, such as reentrant ventricular tachycardia and fibrillation. The mechanisms which lead to the onset of such electrical disorders in patients affected by BrS are not completely understood, yet. The aim of the present study is to investigate by means of numerical simulations the electrophysiological mechanisms at the basis of the morphology of electrocardiogram (ECG) and the onset of reentry associated with BrS. To this end, we consider the Bidomain equations coupled with the ten Tusscher–Panfilov membrane model, on an idealized wedge of human right ventricular tissue. The results have shown that: (1) epicardial dispersion of repolarization, generated by the coexistence of regions of early and late repolarization, due to different modulation of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_\mathrm{{CaL}}$$\end{document}ICaL current, produces ECG waveforms exhibiting qualitatively the typical BrS morphology, characterized by ST elevation and partially negative T-waves; (2) epicardial dispersion of repolarization promotes the onset of reentry during the implementation of the programmed stimulation protocol, because of the conduction block occurring when a premature beat reaches the border of late repolarizing regions; and (3) the modulation of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_\mathrm{{to}}$$\end{document}Ito current affects the duration of reentry, which becomes sustained with a remarkable increase of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_\mathrm{{to}}$$\end{document}Ito in the subepicardial layers. Supplementary Information The online version contains supplementary material available at 10.1007/s11538-023-01124-9.


Introduction
The Brugada syndrome (BrS), first described as a new clinical entity in Brugada and Brugada (1992), has attracted great interest because of its association with the risk of sudden cardiac death and susceptibility to ventricular arrhythmias (VA) in young patients; see Antzelevitch (2001b), Antzelevitch (2006), Benito et al. (2009) andSieira et al. (2016). BrS is an inherited arrhythmogenic disease, transmitted as an autosomal dominant trait, and shows age-and sex-related penetrance. Clinical manifestations of the disease are more frequent in young adults, and they are eightfold more frequent in men than in women. At least 12 genes have been associated with BrS, but only two (SCN5A and CACN1Ac) individually account for more than 5% of positively genotyped patients. Diagnosis of BrS is based on a typical electrocardiographic pattern, characterized by an ST segment elevation of about 2 mm followed by a negative Twave in the right precordial leads (V1-V3), observed either spontaneously or during a sodium channel blocker test (Fowler and Priori 2009;Priori et al. 2002). The incidence of arrhythmic events in patients with BrS was 13.5% per year in patients with a history of sudden cardiac arrest, 3.2% per year in patients with syncope and 1% per year in asymptomatic patients (Fauchier et al. 2013). The mechanisms underlying the characteristic electrocardiographic features and the onset of such arrhythmic events remain poorly understood. For this reason, it is difficult to determine whether a patient is at risk of undergoing arrhythmic events, i.e., to perform the so-called arrhythmic risk stratification. Accurate identification and treatment of individuals at high risk of sudden death are major challenges in the clinical management of BrS patients.
Recent studies in selected patients with Brugada syndrome have described complex arrhythmic substrates in the right ventricular outflow. The development of our work has been inspired by the research published by Nademanee et al. (2017) and Pappone et al. (2018), who explored clinical and electrophysiological predictors of malignant ventricular tachyarrhythmia inducibility in BrS. In particular, they investigated the correlations between the presence and extent of a region of cells with altered electrical properties (malignant VA substrate) in the right ventricle of BrS patients and the pathophysiological basis of lethal VA, which remain still unclear.
Because of the limitation to experimental research involving human cardiac tissue, alternative methods such as computer modeling are of great interest, see Bueno-Orovio et al. (2015), Hoogendijk et al. (2011), Tsumoto et al. (2020) and Xia et al. (2005). Previous computational studies, focusing on the numerical simulation of reentrant ventricular arrhythmias associated with BrS, have considered reduced tissue mod-els, such as the isotropic Monodomain model (Bueno-Orovio et al. 2015), simplified geometries (Tsumoto et al. 2020;Caligari and Scacchi 2020) or phenomenological membrane models (Bueno-Orovio et al. 2015).
The aim of the present work is to investigate the electrophysiological mechanisms at the basis of the morphology of electrocardiogram (ECG) and the onset of reentry associated with BrS, by means of three-dimensional parallel numerical simulations. From the mathematical viewpoint, we have adopted: (1) the anisotropic Bidomain representation of the cardiac tissue, coupled with the human ten Tusscher membrane model; and (2) a parallel finite element solver with GPU acceleration.
The rest of the paper is organized as follows: Sects. 2 and 3 are devoted to the description of the mathematical models and numerical methods adopted for the numerical simulations; Sect. 4 presents the results of the numerical simulations, which are then discussed in Sect. 5.

The Bidomain Model
The macroscopic Bidomain representation of the cardiac tissue volume is obtained by considering the superposition of two anisotropic continuous media, the intra-(i) and extra-(e) cellular media, coexisting at every point of the tissue and separated by a distributed continuous cellular membrane; see, e.g., Neu and Krassowska (1993), Pennacchio et al. (2006) for a derivation of the Bidomain model from homogenization of cellular models and Colli Franzone and Savaré (2002), Bourgault et al. (2009) and Veneroni (2009) for the well-posedness analysis. We recall that the cardiac tissue consists of an arrangement of fibers that rotate counterclockwise from epi-to endocardium, and that have a laminar organization modeled as a set of muscle sheets running radially from epi-to endocardium. The anisotropy of the intra-and extracellular media, related to the macroscopic arrangement of the cardiac myocytes in the fiber structure, is described by the anisotropic conductivity tensors D i (x) and D e (x), respectively, defined in (2).
We denote by ⊂ R 3 the bounded physical region occupied by the cardiac tissue and introduce the parabolic-elliptic formulation of the Bidomain system. Given an applied intracellular current per unit volume I i app : × (0, T ) → R, and initial conditions v 0 : → R, w 0 : → R N w , find the transmembrane potential v : × (0, T ) → R, extracellular potential u e : × (0, T ) → R, the gating and ionic concentrations variables w : where C m is the membrane capacity, χ is the membrane surface to volume ratio and n is the outward unit normal with respect to the domain boundary ∂ . Since the extracellular potential u e is uniquely determined only up to a constant in space, we fix this constant by imposing the condition u e dx = 0. The nonlinear reaction term I ion and the system of ordinary differential equations (ODEs) for the gating and ionic concentration variables w are given by the ionic membrane model; here, we will consider the ten Tusscher (TP06) membrane model (ten Tusscher et al. 2004;ten Tusscher and Panfilov 2006). N w denotes the number of gating and ionic concentration variables, which in the TP06 model amounts to 18.
The conductivity tensors D i (x) and D e (x) at any point x ∈ are defined as Here a l (x), a t (x), a n (x), is a triplet of orthonormal principal axes with a l (x) parallel to the local fiber direction, a t (x) and a n (x) tangent and orthogonal to the radial laminae, respectively, and both being transversal to the fiber axis (see, e.g., LeGrice et al. 1995). Moreover, σ i,e l , σ i,e t , σ i,e n are the conductivity coefficients in the intra-and extracellular media measured along the corresponding directions a l , a t , a n .

Computation of Pseudo-electrocardiograms
In order to compute the extracardiac electrocardiograms (ECGs), we adopt the infinite medium approximation, which consists of assuming that both the bulk medium (intraand extracellular) and extracardiac domains are isotropic media, with the same conductivity coefficient σ b , and the extracardiac domain is unbounded; thus, D i + D e = σ b I . As a consequence, with these simplifications, given a transmembrane potential distribution v inside the cardiac domain , the extracardiac/extracellular potential u in R 3 satisfies the following differential problem: where σ b is the conductivity of the extracardiac/bulk medium. Thus, by exploiting the fundamental solution of the three-dimensional Laplace equation in R 3 , the electric potential u at an extracardiac site x, called pseudo-electrocardiogram (pECG), is computed according to the following formula: For further details, see Colli Franzone et al. (2014, Ch. 5, Proposition 5.1) and Geselowitz and Miller (1983).

Variational Formulation
Let V be the Sobolev space H 1 ( ), define the spaces define the usual L 2 -inner product (ϕ, ψ) = ϕψdx ∀ϕ, ψ ∈ L 2 ( ), and the elliptic bilinear forms where D = D i + D e is the bulk conductivity tensor.

Numerical Methods
In this section, we briefly describe our numerical approximation of the Bidomain model; we refer the interested reader to Vigmond et al. (2002), Colli Franzone and Pavarino (2004), Colli Franzone et al. (2005, 2014, 2018 for further details.

Space Discretization
System (5) is first discretized in space by the finite element method. Let T h be a quasi-uniform triangulation of having maximal diameter h and V h be an associated conforming finite element space. Once a finite element basis {ϕ p } N p=1 of V h is chosen, denote by M = {m pj } the diagonal mass matrix obtained by the usual mass-lumping technique and by A i,e = {a i,e pj } the symmetric intra-and extracellular stiffness matrices, with elements The semi-discrete Bidomain problem, obtained by applying a standard Galerkin procedure, can be written in compact matrix form as with block mass and stiffness matrices ..., w N w ) and I i app , respectively.

Time Discretization
As time discretization, we employ an implicit-explicit (IMEX) strategy, based on decoupling the ODEs from the PDEs and on treating the linear diffusion terms implicitly and the nonlinear reaction terms explicitly. The implicit treatment of the diffusion term is needed in order to avoid a stability constraint on the time step t induced by the fine mesh size h. Nevertheless, due to the explicit treatment of the reaction terms, stability could be preserved for a time step t satisfying a condition of Courant-Friedrichs-Lewy (CFL) type; see, e.g., Quarteroni and Valli (1994). The two equations in (6) arising from the discretization of the PDEs are solved uncoupled. In particular, at the general time step, given w n , v n and u n e computed at the previous time step, -we first update the gating and ionic concentration variables w n+1 by solving the ODEs of the membrane model, -then we solve the elliptic equation computing u n+1 e , -and finally we update the transmembrane potential v n+1 by solving the parabolic equation.
Summarizing in formulae, given w n , v n , u n e , the scheme is As a consequence, at each time step we solve once the linear system with matrix A i + A e arising from the elliptic equation, and once the linear system with matrix c m t M + A i arising from the parabolic equation. Both linear systems are solved by the preconditioned conjugate gradient (PCG) method, since the matrices are symmetric positive definite in the parabolic case and positive semi-definite in the elliptic case. The preconditioner used for the parabolic system is Block Jacobi (BJ), because the related matrix is well conditioned, while for the elliptic system we use an algebraic multigrid preconditioner.

Computational Domain, Discretization and Tissue Parameters
The right ventricular wedge domain is the image of a Cartesian slab using ellipsoidal coordinates, described by the parametric equations 3 cm, c 1 = 5.9 cm, c 2 = 6.4 cm, and φ min = −π/2, φ max = π/2, θ min = −3π/8, θ max = π/8; see Fig. 1. We will refer to the inner surface of the truncated ellipsoid (r = 0) as endocardium and to the outer surface (r = 1) as epicardium. In all computations, a structured grid of 256 × 256 × 24 hexahedral isoparametric Q 1 finite elements of size h ≈ 0.02 cm is used in space, for a total amount of 1 651 225 mesh nodes. Fibers rotate transmurally, linearly with the depth and counterclockwise from epicardium to endocardium, for a total amount of 90 • . For the semi-implicit discretization in time, we use a constant time step size t = 0.05 ms. We also assume a transversely isotropic tissue with conductivity coefficients σ e l = 2, all given in mS. The membrane capacity is set to C m = 1 µF, while the membrane surface to volume ratio is set to χ = 10 3 cm −1 .

BrS Parameter Setting
In the ECG waveforms, a coved ST segment elevation pattern followed by a negative T-wave is considered diagnostic of the Brugada syndrome; see, e.g., Sieira et al. (2016, Fig. 1). A reduced function of cardiac sodium channel plays an important role in the mechanism of BrS: The sodium channel blockers can provoke or augment the Brugada ECG pattern and mutations in SCN5A, the gene encoding the α-subunit of the cardiac sodium channel, are identified in ∼20% of patients (Fowler and Priori 2009;Priori et al. 2002).
The disorder more frequently observed in mutations in SCN5A is a decrease of sodium current I Na (Clancy and Rudy 2002), which leads to an imbalance between the positive inward and outward currents at the end of the transient repolarization phase of the cell action potential. This imbalance should occur with decreases of inward sodium I Na and L-type calcium I CaL currents or an increase of outward potassium I to current, which leads to a development of a characteristic notch and the loss of action potential dome. Thus, one of the ways in which arrhythmias may occur in Brugada syndrome is through the formation of a heterogeneous substrate, in which a region exhibits abnormal action potentials. Such a substrate might be responsible for the onset of reentrant arrhythmias, as suggested in Pappone et al. (2018).
Following the approach of Hoogendijk et al. (2011), where ST segment elevation is connected to reductions of sodium channel (G Na ) and L-type calcium channel (G CaL ) maximal conductances and an increase of potassium channel (G to ) maximal conductance, in our simulations, we modify in the TP06 model G Na , G CaL and G to in the region of BrS cells.
In the following, we will divide the transmural wall of the computational domain into two parts of equal depths, denoted by subendocardial and subepicardial regions, respectively. In order to model transmural heterogeneities of the action potential duration (APD), we scale the I Ks current of a factor 0.7 and 1.4 in the subendocardial and subepicardial regions, respectively. Moreover, we will denote as malignant VA substrate a portion of tissue, located in the basal upper subepicardial region, occupying about 4.5 × 3.5 cm 2 of epicardial surface and 50% transmural depth, as displayed in Fig. 1.
We consider the following different settings in the TP06 model: control: all parameters in the TP06 model are set to their default values as given in the original paper (ten Tusscher and Panfilov 2006); -BRU1: G Na is set in the subepicardial region at 40% of its normal value; -BRU2: in addition to the BRU1 condition, G CaL is set in the subepicardial region at 10% of its normal value; -BRU3: in addition to the BRU2 conditions, G CaL is set in the malignant VA substrate at 120% of its normal value; -BRU4: in addition to the BRU3 conditions, G to is set in the subepicardial region at 400% of its normal value.
The corresponding subepicardial action potential waveforms, computed by solving the ODEs of the TP06 model applying the stimulus at a basic cycle length (BCL) of 500 ms, are reported in Fig. 2.  (BRU3 and BRU4). Note that the waveforms of the BRU2 and BRU3 coincide because the only difference between the two settings is the presence of the malignant VA substrate, thus outside it the cell membrane is the same

Stimulation Protocol
Arrhythmic risk stratification in BrS is still challenging, especially in asymptomatic cases. In general, patients with documented ventricular fibrillation or with unexplained  (Brugada et al. 2003), on the basis of data from their registry, that arrhythmia inducibility at programmed electrical stimulation (PES) could be useful to identify patients at high risk. For our simulations we have decided to use the same pacing procedure introduced in Brugada et al. (2003). Except where otherwise stated, the excitation process is started by applying a stimulus of 350 mA/cm 3 for 1 ms on a small area 0.12×0.12×0.06 cm 3 at the epicardial site indicated in Fig. 1. For each simulation setting, we first apply four pacing stimuli (S1) at a BCL of 500 ms. Then, a premature stimulus (S2) is delivered 350 ms after S1. If S2 does not generate a reentrant arrhythmia, the S1-S2 coupling interval is shortened in steps of 10 ms until arrhythmia is induced or S2 fails to trigger excitation. If arrhythmia is not induced, an additional S3, and if necessary S4 stimulus, is delivered in the same manner as S2 (initially delivered 350 ms after the previous stimulus, and then shortened until arrhythmia is induced or the stimulus fails to induce arrhythmia). The criterion adopted for successful arrhythmia induction is the onset of a reentrant excitation that remains sustained up to the end of simulation time, set to 4 s after the last stimulus. The total cost of a simulation where the reentrant excitation is sustained up to 4 s amounts to approximately 15 h.

Results
We report in this section the results of parallel numerical simulations on the threedimensional wedge of right ventricular tissue. Except where otherwise stated, the simulations have been performed on the Marconi100 and DGX Linux clusters at the CINECA laboratory (see below). Our in-house code is written in C and is based on the parallel libraries MPI and PETSc (Balay et al. 2016), developed at the Argonne National Laboratory (USA). A CUDA kernel has been developed and included in the C code to solve the ODEs of the membrane model on the GPUs. The rest of the section is organized as follows: in Test 1, we study the performance of the parallel solver; in Test 2, we compute the activation, repolarization and APD spatial distribution in the five settings; in Test 3, we study how the different BrS parameter calibrations affect the pECG waveforms; in Test 4, we simulate the induction of reentry.

Test 1: Performance of the GPU and CPU Bidomain Solvers
We compare the performance of the GPU Bidomain solver with our previous CPU solver. Two  In the Marconi100 cluster the GPU computations have been performed using only one core and one GPU, while in the DGX cluster we have used up to eight cores and eight GPUs. On both clusters, the computations have been executed on a single node. In the GPU solver, the parabolic system is preconditioned by the Jacobi preconditioner, while the elliptic system is preconditioned by the PCGAMG algebraic multigrid preconditioner provided by the PETSc library (Balay et al. 2016). In the CPU solver, the parabolic system is preconditioned by the Block Jacobi preconditioner, while the elliptic system is preconditioned by the BoomerAMG algebraic multigrid preconditioner (Henson and Yang 2002) provided by the HYPRE library (https://www.llnl. The results obtained on Marconi100, reported in Table 1, show that: the GPU membrane solver is about 50 times faster than the 64 cores CPU membrane solver; the GPU parabolic solver is about 8 times faster than the 64 cores CPU parabolic solver; the GPU elliptic solver is about twice as fast as the 64 cores CPU elliptic solver; as a result, the GPU Bidomain solver is in total about twice as fast as the 64 cores CPU elliptic solver.
The results obtained on DGX, reported in Table 2, show that: the GPU membrane solver is about 10 times faster than the 64 cores CPU membrane solver; the GPU parabolic solver is about 6 times faster than the 64 cores CPU parabolic solver; the GPU elliptic solver is about 31 times faster than the 64 cores CPU elliptic solver; as a result, the GPU Bidomain solver is in total about 22 times faster than the 64 cores CPU elliptic solver.

Test 2: Activation and Repolarization Sequences and APD Distributions
The aim of this test is to study how the different BrS settings affect the activation and repolarization sequences and the APD distributions, after a unipolar epicardial stimulus applied at the location indicated in Fig. 1. Figures 3, 4, 5, 6 and 7 report the endocardial and epicardial activation time, repolarization time and APD distributions    All quantities are reported in ms in the control, BRU1, BRU2, BRU3 and BRU4 settings, respectively. Table 3 reports the dispersion (maximum minus minimum) of the activation and repolarization time, and APD, computed on endocardium, epicardium and on the entire tissue.
We first observe that, in the BRU1 setting, the activation sequence is about 57 ms longer than in the control setting (239.58 ms (BRU1) vs. 183.99 ms (control)). This delay is due to the reduction of the maximal conductance of the I Na current. The slower activation induces a slower and longer repolarization sequence, that yields, especially on the epicardium, an increase of the APD dispersion of 5 ms (16.6 ms (BRU1) vs. 11.6 ms (control)).
In the BRU2 setting, the activation sequence is comparable with that obtained in the BRU1 setting. Due to the reduction of the maximal conductance of the I CaL current in the subepicardial region, the epicardial APD is about 100 ms shorter than in the BRU1 setting. As a consequence, the epicardial repolarization occurs much earlier than in the BRU1 setting. The electrotonic effect induces a reduction of the repolarization time of about 40 ms on the endocardium, too. However, the endocardial and epicardial dispersions of repolarization remain comparable with those obtained in the BRU1 setting, while the total dispersion of repolarization increases of about 3 ms (298.7 ms (BRU2) vs. 245.5 ms (BRU1)).
In the BRU3 setting, the activation sequence is again comparable with the previous BRU1 and BRU2 settings. The presence of the subepicardial malignant VA substrate, with a larger I CaL current than in the rest of the subepicardial region, yields an increase of the maximal epicardial APD of almost 100 ms with respect to the BRU2 setting (291.22 ms (BRU3) vs. 194.37 ms (BRU2)). As a consequence, the epicardial dispersion of repolarization increases from 219.6 (BRU2) to 250.4 ms (BRU3). Sharp repolarization gradients appear at the boundaries of the malignant VA substrate. The electrotonic effect induces also a modification of the endocardial repolarization sequence, which presents concave isochrones corresponding to the endocardial projection of the malignant VA substrate. Nevertheless, the endocardial and total dispersions of repolarization remain comparable with those obtained in the BRU2 setting.
Finally, the BRU4 setting presents activation, repolarization and APD distribution and dispersions analogous to those of the BRU3 setting.

Test 3: Effects of BrS Model Calibration on the pECG
The aim of this test is to study the role played by the different BrS parameter settings BRU1-4 on the genesis of the ECG morphology. To this end, we compute the pECG according to Eq. (4) in an extracardiac site x located at a distance of 2 cm from the center of the epicardial surface. In these simulations, we mimic the sinus rhythm excitation by stimulating the entire endocardial surface.
In the control setting, the pECG waveform displayed in Fig. 8 exhibits the typical positive QRS complex, since the excitation wavefront propagates from endocardium toward epicardium. Then, after a flat ST interval, the T-wave presents a positive dome, since the subepicardial layers recover before the endocardial ones, due to the transmural APD heterogeneities. In the BRU1 setting, the pECG waveform is comparable with that computed in the control setting, except a longer QRS complex, due to the slower activation induced by the reduction of I Na current.
In the BRU2 setting, the pECG waveform exhibits a large ST elevation and a huge positive T-wave, because of the fast repolarization occurring in the subepicardial layers, where the I CaL current is reduced.
In the BRU3 setting, we observe again a marked ST elevation and an initially positive T-wave, due to the fast repolarization of the apical and central subepicardial layers. However, after about 250 ms, the pECG undergoes an inversion of the T-wave, which becomes negative, because of the late repolarization of the malignant VA substrate. Thus, in this setting, the pECG waveform exhibits the typical characteristics of a BrS ECG, with ST elevation and T-wave inversion.
Finally, in the BRU4 setting, the pECG waveform is analogous to that computed in the BRU3 setting, except a more prominent J-wave. Thus, this result suggests that the increase of I to current in the subepicardial layers seems responsible for the presence of prominent J-waves in BrS ECGs.

Test 4: Induction of Reentry
In this test, we apply the PES protocol with a unipolar stimulating electrode at the epicardial site indicated in Fig. 1. In the control, BRU1 and BRU2 settings, we were not able to induce reentry.
In the BRU3 setting, we applied an S2 stimulus at 260 ms after the S1 stimulus. See movie SM_BRU3 in the Supplementary material. The excitation wavefront elicited from the S2 stimulus propagates through the epicardial surface until it reaches the bottom border of the malignant VA substrate, where a block of excitation occurs. The electric impulse then propagates around the malignant VA substrate, still refractory. When the tissue becomes excitable again, the wavefront enters the malignant VA substrate, generating a reentrant activation, which propagates backward toward the region Fig. 9 Transmembrane and extracellular potential waveforms computed from BRU3 setting after S1-S2 stimulations, in two epicardial sites: one located inside the malignant VA substrate (red line) and one outside (black line), see Fig. 1. t = 0 is the onset of the S1 stimulus. The S2 stimulus is applied at t = 260 ms (Color figure online) of the excitation block and reexcites the epicardial tissue. The reentry is maintained for about 3 s, then it dies; see Fig. 9.
In the BRU4 setting, we applied an S2 stimulus at 260 ms after the S1 stimulus and an S3 stimulus at 220 ms after the S2 stimulus. See Fig. 10 and movie SM_BRU4 in the Supplementary material, starting from the S2 stimulus. The excitation wavefront elicited from the S3 stimulus propagates through the epicardial surface until it reaches the bottom border of the malignant VA substrate, where it is blocked because the tissue is still refractory. When the tissue becomes excitable again, the wavefront enters the malignant VA substrate, generating reentry. The reentry in this case is maintained until the end of the simulation at 4 s; see Fig. 11.

Discussion
By means of parallel three-dimensional simulations, we have studied the electrophysiological mechanisms determining the ECG morphology and the onset of reentry associated with BrS. The cardiac domain considered is a three-dimensional wedge of the right ventricular wall. Despite the simplified geometry, we have taken into account the main features of cardiac electrophysiological modeling, i.e., the anisotropic Bidomain representation of the ventricular tissue, transmural fiber rotation, a mechanistic human ventricular membrane model, and transmural action potential heterogeneities. We have employed the Bidomain model instead of the Monodomain model in order to compute the electrograms inside the cardiac tissue and to have a more accurate representation of the cardiac sources, that might influence significantly the dynamics of reentry in pathological conditions. In order to induce reentry, we have implemented in silico the programmed electrical stimulation protocol, usually adopted in the clinical practice.
The results have shown that: -Epicardial dispersion of repolarization, generated by the coexistence of regions of early and late repolarization (in the malignant VA substrate), due to different modulation of the I CaL current, produces pECG waveforms exhibiting qualitatively the typical BrS morphology, characterized by ST elevation and partially negative T-waves; -Epicardial dispersion of repolarization promotes the onset of reentry during the implementation of the PES protocol, because of the conduction block occurring when a premature beat reaches the border of the late recovering malignant VA substrate; -The 400% increase of I to in the subepicardial layers influences the morphology of the pECG, yielding a prominent J-wave, and the duration of reentry, which becomes sustained.
The presence of subepicardial regions of early repolarization, due to a severe reduction of I CaL current, plays a decisive role in the genesis of ST elevation in the BrS ECG waveform. On the other hand, the concurrent presence of a sufficiently large subepicardial region of late repolarization induces the inversion of T-wave polarity, typical of BrS ECG waveforms, and promotes conduction blocks after premature beats, leading to the onset of reentry. This is in agreement with the experimental investigations (Nademanee et al. 2017;Pappone et al. 2018) and a previous computational study (Bueno-Orovio et al. 2015), where a phenomenological membrane model is adopted, together with a reduced isotropic Monodomain model.
The presence of prominent J-waves in some ECG leads, a characteristic of J-wave syndromes, i.e., a family of syndromes BrS belongs to, has been identified as a marker for a substrate capable of generating life-threatening ventricular arrhythmias, see, e.g., Yan (2010, 2015). Indeed, our computational results are in agreement with these experimental findings, because we have shown that large I to values yield prominent J-waves and sustained reentry.

Clinical Implications
The results of this study confirm the rationale of selective pharmacological modulation of I to and I CaL currents and/or of transcatheter ablation of the right ventricular epicardial substrate to reduce the likelihood of ventricular arrhythmias of BrS patients.

Limitations
A simplified geometric wedge model of the right ventricular tissue has been employed in this study. The extension to a biventricular geometry, with the inclusion of a patientspecific characterization of the malignant substrate based on electro-anatomical mapping, would strengthen the results obtained in the present study. A further limitation of the study is the computation of the pECG instead of coupling the Bidomain equations with the Laplace equation in the surrounding medium. The rationale for this choice was the need of reducing the computational costs. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.